library(plyr)
library(igraph)
library(RgoogleMaps)
library(RColorBrewer)
library(fmsb)
library(Hmisc)
library(MASS)
library(stargazer)
library(tidyverse)
library(multcomp)
library(WDI)
library(countrycode)

setwd("~/Dropbox/China Huawei Paper/Production Materials/Replication/")
master.data <- read.csv("dataFile_China_Huawei.csv"); head(master.data)


# 250,000
dat <- master.data
dat$treated = 0
threshold=250000

dat = as.data.frame(dat %>%
  group_by(cabb) %>%
  dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
         treated = ifelse(china.huawei > threshold, 1, treated)) %>%
  fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0

dat$ptBOL = 1*(dat$cabb=="BOL")*dat$treated
dat$ptCDI = 1*(dat$cabb=="CDI")*dat$treated
dat$ptCOS = 1*(dat$cabb=="COS")*dat$treated
dat$ptGHA = 1*(dat$cabb=="GHA")*dat$treated
dat$ptINS = 1*(dat$cabb=="INS")*dat$treated
dat$ptMAW = 1*(dat$cabb=="MAW")*dat$treated
dat$ptMLD = 1*(dat$cabb=="MLD")*dat$treated
dat$ptNIG = 1*(dat$cabb=="NIG")*dat$treated
dat$ptUKR = 1*(dat$cabb=="UKR")*dat$treated
dat$ptURU = 1*(dat$cabb=="URU")*dat$treated
dat$ptVEN = 1*(dat$cabb=="VEN")*dat$treated

dat = dat[dat$v2x_polyarchy>0.42,]

m1 <- lm(v2smgovfilprc ~ ptBOL + ptCDI + ptCOS + ptGHA + ptINS + ptMAW + ptMLD + ptNIG + ptUKR + ptURU + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m1)

m2 <- lm(v2smgovshut ~ ptBOL + ptCDI + ptCOS + ptGHA + ptINS + ptMAW + ptMLD + ptNIG + ptUKR + ptURU + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m2)

m3 <- lm(v2smgovsmmon ~ ptBOL + ptCDI + ptCOS + ptGHA + ptINS + ptMAW + ptMLD + ptNIG + ptUKR + ptURU + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m3)

m4 <- lm(v2smarrest ~ ptBOL + ptCDI + ptCOS + ptGHA + ptINS + ptMAW + ptMLD + ptNIG + ptUKR + ptURU + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m4)

n.treated <- length(sort(unique(dat$cabb[dat$treated==1]))) - 4
cont <- matrix(c(0, sum(dat$ptBOL,na.rm=T), sum(dat$ptCDI,na.rm=T), sum(dat$ptCOS,na.rm=T), sum(dat$ptGHA,na.rm=T), sum(dat$ptINS,na.rm=T), sum(dat$ptMAW,na.rm=T), sum(dat$ptMLD,na.rm=T), sum(dat$ptNIG,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptURU,na.rm=T), sum(dat$ptVEN,na.rm=T), rep(0, (length(coef(m1))-n.treated-1) )),1)

summary(glht(m1,cont))
summary(glht(m2,cont))
summary(glht(m3,cont))
summary(glht(m4,cont))

total.treatment.periods <- sum(sum(dat$ptBOL,na.rm=T), sum(dat$ptCDI,na.rm=T), sum(dat$ptCOS,na.rm=T), sum(dat$ptGHA,na.rm=T), sum(dat$ptINS,na.rm=T), sum(dat$ptMAW,na.rm=T), sum(dat$ptMLD,na.rm=T), sum(dat$ptNIG,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptURU,na.rm=T), sum(dat$ptVEN,na.rm=T))
atreat.m1 <- glht(m1,cont/total.treatment.periods); summary(atreat.m1)
atreat.m2 <- glht(m2,cont/total.treatment.periods); summary(atreat.m2)
atreat.m3 <- glht(m3,cont/total.treatment.periods); summary(atreat.m3)
atreat.m4 <- glht(m4,cont/total.treatment.periods); summary(atreat.m4)


m1.250coef = summary(atreat.m1)$test[3]$coefficients 
m1.250se = summary(atreat.m1)$test[4]$sigma 
m1.250p = summary(atreat.m1)$test[6]$pvalues 

m2.250coef = summary(atreat.m2)$test[3]$coefficients 
m2.250se = summary(atreat.m2)$test[4]$sigma 
m2.250p = summary(atreat.m2)$test[6]$pvalues 

m3.250coef = summary(atreat.m3)$test[3]$coefficients 
m3.250se = summary(atreat.m3)$test[4]$sigma 
m3.250p = summary(atreat.m3)$test[6]$pvalues 

m4.250coef = summary(atreat.m4)$test[3]$coefficients 
m4.250se = summary(atreat.m4)$test[4]$sigma 
m4.250p = summary(atreat.m4)$test[6]$pvalues 

# 500,000
dat <- master.data
dat$treated = 0
threshold=500000

dat = as.data.frame(dat %>%
  group_by(cabb) %>%
  dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
         treated = ifelse(china.huawei > threshold, 1, treated)) %>%
  fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0

dat$ptBOL = 1*(dat$cabb=="BOL")*dat$treated
dat$ptCDI = 1*(dat$cabb=="CDI")*dat$treated
dat$ptCOS = 1*(dat$cabb=="COS")*dat$treated
dat$ptGHA = 1*(dat$cabb=="GHA")*dat$treated
dat$ptINS = 1*(dat$cabb=="INS")*dat$treated
dat$ptMAW = 1*(dat$cabb=="MAW")*dat$treated
dat$ptMLD = 1*(dat$cabb=="MLD")*dat$treated
dat$ptUKR = 1*(dat$cabb=="UKR")*dat$treated
dat$ptURU = 1*(dat$cabb=="URU")*dat$treated
dat$ptVEN = 1*(dat$cabb=="VEN")*dat$treated

dat = dat[dat$v2x_polyarchy>0.42,]

m1 <- lm(v2smgovfilprc ~ ptBOL + ptCDI + ptCOS + ptGHA + ptINS + ptMAW + ptMLD + ptUKR + ptURU + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m1)

m2 <- lm(v2smgovshut ~ ptBOL + ptCDI + ptCOS + ptGHA + ptINS + ptMAW + ptMLD + ptUKR + ptURU + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m2)

m3 <- lm(v2smgovsmmon ~ ptBOL + ptCDI + ptCOS + ptGHA + ptINS + ptMAW + ptMLD + ptUKR + ptURU + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m3)

m4 <- lm(v2smarrest ~ ptBOL + ptCDI + ptCOS + ptGHA + ptINS + ptMAW + ptMLD + ptUKR + ptURU + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m4)

n.treated <- length(sort(unique(dat$cabb[dat$treated==1]))) - 4
cont <- matrix(c(0, sum(dat$ptBOL,na.rm=T), sum(dat$ptCDI,na.rm=T), sum(dat$ptCOS,na.rm=T), sum(dat$ptGHA,na.rm=T), sum(dat$ptINS,na.rm=T), sum(dat$ptMAW,na.rm=T), sum(dat$ptMLD,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptURU,na.rm=T), sum(dat$ptVEN,na.rm=T), rep(0, (length(coef(m1))-n.treated-1) )),1)

summary(glht(m1,cont))
summary(glht(m2,cont))
summary(glht(m3,cont))
summary(glht(m4,cont))

total.treatment.periods <- sum(sum(dat$ptBOL,na.rm=T), sum(dat$ptCDI,na.rm=T), sum(dat$ptCOS,na.rm=T), sum(dat$ptGHA,na.rm=T), sum(dat$ptINS,na.rm=T), sum(dat$ptMAW,na.rm=T), sum(dat$ptMLD,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptURU,na.rm=T), sum(dat$ptVEN,na.rm=T))
atreat.m1 <- glht(m1,cont/total.treatment.periods); summary(atreat.m1)
atreat.m2 <- glht(m2,cont/total.treatment.periods); summary(atreat.m2)
atreat.m3 <- glht(m3,cont/total.treatment.periods); summary(atreat.m3)
atreat.m4 <- glht(m4,cont/total.treatment.periods); summary(atreat.m4)

m1.500coef = summary(atreat.m1)$test[3]$coefficients 
m1.500se = summary(atreat.m1)$test[4]$sigma 
m1.500p = summary(atreat.m1)$test[6]$pvalues 

m2.500coef = summary(atreat.m2)$test[3]$coefficients 
m2.500se = summary(atreat.m2)$test[4]$sigma 
m2.500p = summary(atreat.m2)$test[6]$pvalues 

m3.500coef = summary(atreat.m3)$test[3]$coefficients 
m3.500se = summary(atreat.m3)$test[4]$sigma 
m3.500p = summary(atreat.m3)$test[6]$pvalues 

m4.500coef = summary(atreat.m4)$test[3]$coefficients 
m4.500se = summary(atreat.m4)$test[4]$sigma 
m4.500p = summary(atreat.m4)$test[6]$pvalues 


# 1 million
dat <- master.data
dat$treated = 0
threshold=1000000

dat = as.data.frame(dat %>%
  group_by(cabb) %>%
  dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
         treated = ifelse(china.huawei > threshold, 1, treated)) %>%
  fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0

dat$ptBOL = 1*(dat$cabb=="BOL")*dat$treated
dat$ptCDI = 1*(dat$cabb=="CDI")*dat$treated
dat$ptCOS = 1*(dat$cabb=="COS")*dat$treated
dat$ptGHA = 1*(dat$cabb=="GHA")*dat$treated
dat$ptINS = 1*(dat$cabb=="INS")*dat$treated
dat$ptMAW = 1*(dat$cabb=="MAW")*dat$treated
dat$ptMLD = 1*(dat$cabb=="MLD")*dat$treated
dat$ptUKR = 1*(dat$cabb=="UKR")*dat$treated
dat$ptVEN = 1*(dat$cabb=="VEN")*dat$treated

dat = dat[dat$v2x_polyarchy>0.42,]

m1 <- lm(v2smgovfilprc ~ ptBOL + ptCDI + ptCOS + ptGHA + ptINS + ptMAW + ptMLD + ptUKR + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m1)

m2 <- lm(v2smgovshut ~ ptBOL + ptCDI + ptCOS + ptGHA + ptINS + ptMAW + ptMLD + ptUKR + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m2)

m3 <- lm(v2smgovsmmon ~ ptBOL + ptCDI + ptCOS + ptGHA + ptINS + ptMAW + ptMLD + ptUKR + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m3)

m4 <- lm(v2smarrest ~ ptBOL + ptCDI + ptCOS + ptGHA + ptINS + ptMAW + ptMLD + ptUKR + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m4)

n.treated <- length(sort(unique(dat$cabb[dat$treated==1]))) - 3
cont <- matrix(c(0, sum(dat$ptBOL,na.rm=T), sum(dat$ptCDI,na.rm=T), sum(dat$ptCOS,na.rm=T), sum(dat$ptGHA,na.rm=T), sum(dat$ptINS,na.rm=T), sum(dat$ptMAW,na.rm=T), sum(dat$ptMLD,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptVEN,na.rm=T), rep(0, (length(coef(m1))-n.treated-1) )),1)

summary(glht(m1,cont))
summary(glht(m2,cont))
summary(glht(m3,cont))
summary(glht(m4,cont))

total.treatment.periods <- sum(sum(dat$ptBOL,na.rm=T), sum(dat$ptCDI,na.rm=T), sum(dat$ptCOS,na.rm=T), sum(dat$ptGHA,na.rm=T), sum(dat$ptINS,na.rm=T), sum(dat$ptMAW,na.rm=T), sum(dat$ptMLD,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptVEN,na.rm=T))
atreat.m1 <- glht(m1,cont/total.treatment.periods); summary(atreat.m1)
atreat.m2 <- glht(m2,cont/total.treatment.periods); summary(atreat.m2)
atreat.m3 <- glht(m3,cont/total.treatment.periods); summary(atreat.m3)
atreat.m4 <- glht(m4,cont/total.treatment.periods); summary(atreat.m4)

m1.1mcoef = summary(atreat.m1)$test[3]$coefficients 
m1.1mse = summary(atreat.m1)$test[4]$sigma 
m1.1mp = summary(atreat.m1)$test[6]$pvalues 

m2.1mcoef = summary(atreat.m2)$test[3]$coefficients 
m2.1mse = summary(atreat.m2)$test[4]$sigma 
m2.1mp = summary(atreat.m2)$test[6]$pvalues 

m3.1mcoef = summary(atreat.m3)$test[3]$coefficients 
m3.1mse = summary(atreat.m3)$test[4]$sigma 
m3.1mp = summary(atreat.m3)$test[6]$pvalues 

m4.1mcoef = summary(atreat.m4)$test[3]$coefficients 
m4.1mse = summary(atreat.m4)$test[4]$sigma 
m4.1mp = summary(atreat.m4)$test[6]$pvalues 


# 5 million
dat <- master.data
dat$treated = 0
threshold=5000000

dat = as.data.frame(dat %>%
  group_by(cabb) %>%
  dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
         treated = ifelse(china.huawei > threshold, 1, treated)) %>%
  fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0

dat$ptBOL = 1*(dat$cabb=="BOL")*dat$treated
dat$ptCDI = 1*(dat$cabb=="CDI")*dat$treated
dat$ptCOS = 1*(dat$cabb=="COS")*dat$treated
dat$ptINS = 1*(dat$cabb=="INS")*dat$treated
dat$ptMAW = 1*(dat$cabb=="MAW")*dat$treated
dat$ptMLD = 1*(dat$cabb=="MLD")*dat$treated
dat$ptUKR = 1*(dat$cabb=="UKR")*dat$treated
dat$ptVEN = 1*(dat$cabb=="VEN")*dat$treated

dat = dat[dat$v2x_polyarchy>0.42,]

m1 <- lm(v2smgovfilprc ~ ptBOL + ptCDI + ptCOS + ptINS + ptMAW + ptMLD + ptUKR + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m1)

m2 <- lm(v2smgovshut ~ ptBOL + ptCDI + ptCOS + ptINS + ptMAW + ptMLD + ptUKR + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m2)

m3 <- lm(v2smgovsmmon ~ ptBOL + ptCDI + ptCOS + ptINS + ptMAW + ptMLD + ptUKR + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m3)

m4 <- lm(v2smarrest ~ ptBOL + ptCDI + ptCOS + ptINS + ptMAW + ptMLD + ptUKR + ptVEN + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m4)

n.treated <- length(sort(unique(dat$cabb[dat$treated==1]))) - 3
cont <- matrix(c(0, sum(dat$ptBOL,na.rm=T), sum(dat$ptCDI,na.rm=T), sum(dat$ptCOS,na.rm=T), sum(dat$ptINS,na.rm=T), sum(dat$ptMAW,na.rm=T), sum(dat$ptMLD,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptVEN,na.rm=T), rep(0, (length(coef(m1))-n.treated-1) )),1)

summary(glht(m1,cont))
summary(glht(m2,cont))
summary(glht(m3,cont))
summary(glht(m4,cont))

total.treatment.periods <- sum(sum(dat$ptBOL,na.rm=T), sum(dat$ptCDI,na.rm=T), sum(dat$ptCOS,na.rm=T), sum(dat$ptINS,na.rm=T), sum(dat$ptMAW,na.rm=T), sum(dat$ptMLD,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptVEN,na.rm=T))
atreat.m1 <- glht(m1,cont/total.treatment.periods); summary(atreat.m1)
atreat.m2 <- glht(m2,cont/total.treatment.periods); summary(atreat.m2)
atreat.m3 <- glht(m3,cont/total.treatment.periods); summary(atreat.m3)
atreat.m4 <- glht(m4,cont/total.treatment.periods); summary(atreat.m4)

m1.5mcoef = summary(atreat.m1)$test[3]$coefficients 
m1.5mse = summary(atreat.m1)$test[4]$sigma 
m1.5mp = summary(atreat.m1)$test[6]$pvalues 

m2.5mcoef = summary(atreat.m2)$test[3]$coefficients 
m2.5mse = summary(atreat.m2)$test[4]$sigma 
m2.5mp = summary(atreat.m2)$test[6]$pvalues 

m3.5mcoef = summary(atreat.m3)$test[3]$coefficients 
m3.5mse = summary(atreat.m3)$test[4]$sigma 
m3.5mp = summary(atreat.m3)$test[6]$pvalues 

m4.5mcoef = summary(atreat.m4)$test[3]$coefficients 
m4.5mse = summary(atreat.m4)$test[4]$sigma 
m4.5mp = summary(atreat.m4)$test[6]$pvalues 



# 10 million
dat <- master.data
dat$treated = 0
threshold=10000000

dat = as.data.frame(dat %>%
  group_by(cabb) %>%
  dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
         treated = ifelse(china.huawei > threshold, 1, treated)) %>%
  fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0

dat$ptBOL = 1*(dat$cabb=="BOL")*dat$treated
dat$ptCDI = 1*(dat$cabb=="CDI")*dat$treated
dat$ptCOS = 1*(dat$cabb=="COS")*dat$treated
dat$ptINS = 1*(dat$cabb=="INS")*dat$treated
dat$ptMAW = 1*(dat$cabb=="MAW")*dat$treated
dat$ptMLD = 1*(dat$cabb=="MLD")*dat$treated
dat$ptUKR = 1*(dat$cabb=="UKR")*dat$treated

dat = dat[dat$v2x_polyarchy>0.42,]

m1 <- lm(v2smgovfilprc ~ ptBOL + ptCDI + ptCOS + ptINS + ptMAW + ptMLD + ptUKR + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m1)

m2 <- lm(v2smgovshut ~ ptBOL + ptCDI + ptCOS + ptINS + ptMAW + ptMLD + ptUKR  + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m2)

m3 <- lm(v2smgovsmmon ~ ptBOL + ptCDI + ptCOS + ptINS + ptMAW + ptMLD + ptUKR + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m3)

m4 <- lm(v2smarrest ~ ptBOL + ptCDI + ptCOS + ptINS + ptMAW + ptMLD + ptUKR + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m4)

n.treated <- length(sort(unique(dat$cabb[dat$treated==1]))) - 1
cont <- matrix(c(0, sum(dat$ptBOL,na.rm=T), sum(dat$ptCDI,na.rm=T), sum(dat$ptCOS,na.rm=T), sum(dat$ptINS,na.rm=T), sum(dat$ptMAW,na.rm=T), sum(dat$ptMLD,na.rm=T), sum(dat$ptUKR,na.rm=T), rep(0, (length(coef(m1))-n.treated-1) )),1)

summary(glht(m1,cont))	#yes
summary(glht(m2,cont))	#yes
summary(glht(m3,cont))	#no
summary(glht(m4,cont))	#yes

total.treatment.periods <- sum(sum(dat$ptBOL,na.rm=T), sum(dat$ptCDI,na.rm=T), sum(dat$ptCOS,na.rm=T), sum(dat$ptINS,na.rm=T), sum(dat$ptMAW,na.rm=T), sum(dat$ptMLD,na.rm=T), sum(dat$ptUKR,na.rm=T))
atreat.m1 <- glht(m1,cont/total.treatment.periods); summary(atreat.m1)
atreat.m2 <- glht(m2,cont/total.treatment.periods); summary(atreat.m2)
atreat.m3 <- glht(m3,cont/total.treatment.periods); summary(atreat.m3)
atreat.m4 <- glht(m4,cont/total.treatment.periods); summary(atreat.m4)


m1.10mcoef = summary(atreat.m1)$test[3]$coefficients 
m1.10mse = summary(atreat.m1)$test[4]$sigma 
m1.10mp = summary(atreat.m1)$test[6]$pvalues 

m2.10mcoef = summary(atreat.m2)$test[3]$coefficients 
m2.10mse = summary(atreat.m2)$test[4]$sigma 
m2.10mp = summary(atreat.m2)$test[6]$pvalues 

m3.10mcoef = summary(atreat.m3)$test[3]$coefficients 
m3.10mse = summary(atreat.m3)$test[4]$sigma 
m3.10mp = summary(atreat.m3)$test[6]$pvalues 

m4.10mcoef = summary(atreat.m4)$test[3]$coefficients 
m4.10mse = summary(atreat.m4)$test[4]$sigma 
m4.10mp = summary(atreat.m4)$test[6]$pvalues 


# Export table

ti = data.frame(rbind(c(m1.250coef, m1.250se, m1.250p),c(m1.500coef, m1.500se, m1.500p),c(m1.1mcoef, m1.1mse, m1.1mp),c(m1.5mcoef, m1.5mse, m1.5mp),c(m1.10mcoef, m1.10mse, m1.10mp)))
ti$term = c("Transfer \u2265 $250,000","Transfer \u2265 $500,000","Transfer \u2265 $1,000,000","Transfer \u2265 $5,000,000","Transfer \u2265 $10,000,000")
colnames(ti) = c("estimate","std.error","p.value","term")
ti = ti[c("term","estimate","std.error","p.value")]
out1 = list(tidy=ti)
class(out1) = "modelsummary_list"

ti = data.frame(rbind(c(m2.250coef, m2.250se, m2.250p), c(m2.500coef, m2.500se, m2.500p), c(m2.1mcoef, m2.1mse, m2.1mp), c(m2.5mcoef, m2.5mse, m2.5mp),c(m2.10mcoef, m2.10mse, m2.10mp)))
ti$term = c("Transfer \u2265 $250,000","Transfer \u2265 $500,000","Transfer \u2265 $1,000,000","Transfer \u2265 $5,000,000","Transfer \u2265 $10,000,000")
colnames(ti) = c("estimate","std.error","p.value","term")
ti = ti[c("term","estimate","std.error","p.value")]
out2 = list(tidy=ti)
class(out2) = "modelsummary_list"


ti = data.frame(rbind(c(m3.250coef, m3.250se, m3.250p),c(m3.500coef, m3.500se, m3.500p),c(m3.1mcoef, m3.1mse, m3.1mp),c(m3.5mcoef, m3.5mse, m3.5mp),c(m3.10mcoef, m3.10mse, m3.10mp)))
ti$term = c("Transfer \u2265 $250,000","Transfer \u2265 $500,000","Transfer \u2265 $1,000,000","Transfer \u2265 $5,000,000","Transfer \u2265 $10,000,000")
colnames(ti) = c("estimate","std.error","p.value","term")
ti = ti[c("term","estimate","std.error","p.value")]
out3 = list(tidy=ti)
class(out3) = "modelsummary_list"

ti = data.frame(rbind(c(m4.250coef, m4.250se, m4.250p),c(m4.500coef, m4.500se, m4.500p),c(m4.1mcoef, m4.1mse, m4.1mp),c(m4.5mcoef, m4.5mse, m4.5mp),c(m4.10mcoef, m4.10mse, m4.10mp)))
ti$term = c("Transfer \u2265 $250,000","Transfer \u2265 $500,000","Transfer \u2265 $1,000,000","Transfer \u2265 $5,000,000","Transfer \u2265 $10,000,000")
colnames(ti) = c("estimate","std.error","p.value","term")
ti = ti[c("term","estimate","std.error","p.value")]
out4 = list(tidy=ti)
class(out4) = "modelsummary_list"


myrows = data.frame(list(c("Country Fixed Effects",rep("Yes",4))
                         ,c("Year Fixed Effects",rep("Yes",4))
                         ,c("Control Variables",rep("Yes",4))
))

export = modelsummary(list("Internet Filtering" = out1
                           , "Internet Shutdowns" = out2
                           ,"Social Media Monitoring"=out3
                           ,"Arrests for Political Content"=out4
)
,add_rows=data.frame(t(myrows))
, stars=c('*' = .1, '**' = 0.05, '***' = .01)
, output="latex_tabular"
)
export = gsub("≥", "$\\\\ge$", export)
export

